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1 Introduction and motivation 

We present a boundary integral method, and an accompanying boundary 
element discretization, for solving boundary-value problems for the Laplace- 
Beltrami operator on the surface of the unit sphere <S in M 3 . We consider a 
closed curve C on S which divides S into two parts Si and 1S2. In particular, 
C = dSi is the boundary curve of Si. We are interested in solving a boundary 
value problem for the Laplace-Beltrami operator in S2, with boundary data 
prescribed on C. 

We shall begin by describing a physical problem of interest. Then, we derive 
an integral representation formula for solutions of the Laplace-Beltrami op- 
erator on the sphere, and introduce the single and double layer potentials. 
We investigate their jump properties, and use these to derive an integral 
equation for the solution of a Dirichlet problem. A variational strategy is 
presented, along with some numerical experiments validating our ideas. To 
the best of our knowledge, these integral equations and their discretizations 
have not been studied before. We believe they present an elegant and natural 
solution strategy for boundary value problems on the sphere. 
This work is motivated in part by recent investigations into the motion 
of point vortices on spheres, specifically in bounded regions with walls on 
the sphere. Kidambi and Newton [4J considered such a problem, assuming 
the bounded sub-surface of the sphere lent itself to the method of images. 
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Crowdy, in a series of papers [H [21 [3], has also investigated the motion 
of vortices on spheres. In [3], he uses conformal mapping onto the complex 
plane to study the motion of a vortex on a sphere with walls. We shall study 
a closely related model problem, for which the methods of O d] would be 
applicable. However, we propose an integral-equation method instead which 
is valid for any bounded sub-region 52, provided the curve C is sufficiently 
smooth. This technique will be valid even where the method of images is 
not, and which does not involve explicit knowledge of conformal mappings 
between the stereographically-projected subregion of interest, and the upper 
half of the complex plane. 

1.1 Point vortex motion on a sphere with walls 

The underlying physical phenomenon considered in [3] is the motion of a 
point vortex in an incompressible fluid on the surface of the unit sphere, 5. 
There is a bounded solid region, denoted <Si C 5, with a simply connected 
boundary, C. No fluid can penetrate into S±. Let 52 be the surface of the 
sphere excluding 5i U C, see Figure (1). 



Figure 1: The unit sphere, 5, with an impenetrable island 5i on the surface. 
C is the boundary of the island. 

A point on the sphere 5 will be described in terms of the spherical angles, 
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We consider a point vortex of strength k located at a point xo G £2- The 
flow motion is assumed irrotational, except for the point vorticity associated 
with the vortex. This assumption needs some justification, which we will 
discuss below. The incompressible nature of the fluid allows us to prescribe 
a stream function, ^(xq, x), for the fluid velocity. That is, the velocity field 
satisfies 

u = x e r . 

Here e r is the unit radial vector to the surface. The vorticity is then defined 

as 

uj = uje r := V x u. 

If the fluid motion is irrotational except at xo, then ui = except at that 
point. We insist that the boundary C be a streamline of the motion. Without 
loss of generality, we can set the streamline constant to zero. The function 
^ is really the Green's function for the Laplace-Beltrami operator on the 
subsurface 52 of the sphere: 

— As ^(xq, x) = k5(\x — xo\), \/x G S2, (1.1a) 

*(xo,x) = 0, ViGC. (1.1b) 

We recall that, in spherical coordinates, A5 is defined as 

1 d 2 1 d ( . „ d 



A s u(x) 



_|_ ( CJITI 

_sin 2 0cV sin 00\ " 39 



u(x(<p,0)). (1.2) 



The assumption that the fluid motion is irrotational can be justified by 
noting that since the fluid is incompressible, it can perfectly slip at S. This 
allows us to prescribe a circulation at C, so that the Gauss constraint for 
the vorticity, 

/ uj da = 0, 
Js 2 

is satisfied. 

Analogously to finding Green's functions in the plane, we can find ^ in 
terms of the fundamental singularity U of the Laplace-Beltrami operator on 
the entire surface of the sphere, and a smooth function v xo (x). That is, 

V(x,x ) = U(x,x ) + v X0 (x) (1.3) 

where v xo solves 

As v xo (x) = VxGfi, v xo (x) = -U(xq,x), xeC. (1.4) 
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We can interpret the fundamental singularity, U, as the stream function for 
a point vortex of strength k on the sphere without boundaries. We denote 
the fundamental singularity, with k = 1, as U* henceforth. Note that from 
[3], we get that 



U(x a , x) = — ft log 



(z - z a )(z - z a ) 



;i + |z| 2 )(i + |z a | 2 ) 



where we've stated this in terms of the mapped points z, z a in the complex 
plane: 

z = cot(0/2)e^, x = {6,<f>) 
The fundamental singularity U satisfies the partial differential equation: 

- A S U = k (d(\x - x \) - ^) , (1.5a) 

and the Gauss condition for the vorticity, u = — Ag U: 

Lods = 0. (1.5b) 
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The implication of (|1.5ap is that there is a "sea" of uniform vorticity, j-, 
in which the point vortex at xq must be embedded if moving on the whole 
sphere. We cannot find a distribution U which satisfies — As U = 5\x — xq\ 
on the sphere, and which simultaneously satisfies the Gauss constraint. In 
order to satisfy the constraint and simultaneously have an irrotational flow, 
we must either counterbalance the point vortex at xq by another vortex 
on the sphere, or have the entire fluid moving with a uniform background 
vorticity. This feature of the fundamental singularity will appear again in 
the next section, and will require us to impose a side constraint on the 
solution density, when employing integral equations. This is reminiscent of 
similar problems arising in the solution of potential problems in unbounded 
regions of the plane. 

At this juncture, we could use the Green's function ^ to study the motion 
of the point vortex, which is governed in the stereographic coordinates z by 



— - T (i + N )^H^ - 



The solutions of this evolution equation are described in terms of level sets 
of the smooth part, v xo (xq) = constant. Such an investigation is performed 
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in the papers by Kidambi and Newton [4j, and Crowdy [3]. Instead, we shall 
study a closely related mathematical model problem. 

Consider the Dirichlet boundary value problem in 52 for the Laplace-Beltrami 
operator: 

Find a smooth u such that for given Dirichlet data g 



We wish to solve (jl.6p by reformulating the boundary value problem as an 
integral equation. As usual, the process of reformulation is not unique; we 
shall be employing a layer ansatz, and solving an integral equation of the 
first kind for the unknown density. We note that we could equivalently have 
chosen to study the Neumann or Robin problem for the system. We could 
use the Green's function for 52, ^f, to solve this Dirichlet problem. The 
methods suggested in [1] and [5] would also be applicable for our model 
problem, with some caveats: the method of Kidambi and Newton relies 
on the ability to use the method of images, while Crowdy's work requires 
knowledge of a conformal map from the stereographically-projected 52 into 
the upper half plane or the unit circle. Instead, we propose an integral 
equation method which is valid for any bounded sub-region 52, provided 
the curve C is sufficiently smooth, and without conformally mapping to the 
plane. Additionally, integral equations allow us to solve problems with lower 
regularity properties, a feature we shall explore in upcoming work. 
If 5i were to degenerate, ie, if the interior of C had zero area, we would need 
to add extra conditions to satisfy the Gauss constraint. Mathematically, we 
would be dealing with the screen problem, and anticipate singular behaviour 
on the corners of the screen. On the surface of the entire sphere without 
walls, we must either embed the point vortex in a fluid of uniform vorticity 
(hence no longer irrotational) , or counter-balance it by another point vortex. 

2 An integral representation formula on the sphere 

We begin by reminding the reader of some vectorial identities on the sphere. 
Let e r ,eQ,e v be the usual unit vectors in spherical coordinates. Recall that 
we can define the surface gradient of a scalar / on 5 as 





(1.6a) 
(1.6b) 




for x £ C 



Vsf(x) 



1 df ^ + df 
sin 9 dip v 89 
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In the same way we introduce the surface divergence for a vector-valued 
function V on the sphere as 

We easily see the identity: 

Agu(x) = div l 5V l su(x). 
We introduce the vectorial surface rotation for a scalar field / on the sphere: 

89 sin 9 8ip 
and the (scalar) surface rotation of a vector field V as 

calls9 ^ = ES (-£ v '<*»> + |o»«Wfc«)>) ■ 

We then obtain another vectorial identity for the Laplace-Beltrami operator: 

Agu(x) = — cutIs cutI s u(x) for x £ S. 

We shall be using a variational setting for most of this paper; to this end, 
we introduce the inner product 

(u,v) L2 ( S ) = J u(x)v(x)da x = J J u(x(ip, 9))v(x(p, 9)) sin 9 d9dp. 

S 

We shall now derive the Green's identiy. We find, by integration by parts, 

{-A s u,v) L2{S) = a s (u,v) = a s {v,u) = (u,-A s v) L ^ S ) 
where we have introduced the symmetric bilinear form 

2tt TV 

1 r) f) f) f) 

d9dip 



a s (u,v) := J J 
o o 



1 8 8 8 8 

u(ip, 9)—v(ip, 9) + sin0— v(ip, 9)—u(ip, 9) 



sin 9 dip dp 89 89 

j Vsu(x) ■ Vsv(x) da. 



6 



Stoke's theorem for the positively oriented curve C and region S2 may be 
written as 

J cwc\sV(x)dcr x = J V(x) • t[x) ds x . 

S 2 C 

Here, t is the unit tangent vector to C. We note that a similar identity holds 
for the region Si, with care taken with the orientation of the tangent. Now, 
setting V = v(x)W(x) and applying the product rule we get 

J curl s v(x) -W(x) da x = — J v(x)[W(x) ■t(x)]ds x + J v(x)cvx\sW(x)da x . 
S2 C 1S2 

With W(x) = curl y it(x) we finally obtain Green's first formula for the 
Laplace-Beltrami operator, 

cutI ^v(x) ■ cui\ s u(x) da x = J v(x)[cm\ s u(x) ■ t(x)]ds x 
s 2 c 



+ J v(x)A s u(x)da x . (2.1) 



5_. 

Note that the left hand side of (|2.ip coincides with the bilinear form as(u, v), 
with the role of S being played by £2 • 

2.1 Fundamental solution and a representation formula 

Proposition 2.1 [HI [7] The fundamental solution of the Laplace-Beltrami 
operator As as defined in (|1 .2j) is given by 

U*(x,x ) = ~log|l-(x,x )| (2.2) 

= — — log [1 — cos((p — tpo) sin 8 sin 6q — cos 6 cos 6q] . (2-3) 
In particular, 

A s U*(x),x ) = 1- (2.4) 
for x = x(ip,8), x = x((p ,8 ) G S with (cp,6) ^ (<po,0 ). 
Remark 2.1 For x,xq € S we have 

\x - x \ 2 = \x\ 2 - 2(x,x ) + \x \ 2 = 2[1 - (x,Xq)]. 
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Hence we obtain 

-— \og\x - x \ = —T- [log[l - (x,x )] + log2] . 

In particular, the fundamental solution of the three-dimensional Laplace- 
Beltrami operator corresponds to the fundamental solution of the two-dimen- 
sional Laplace operator. 

The first Green's identity can be used to derive a representation formula for 
smooth functions defined on the sphere. 

Proposition 2.2 Every sufficiently smooth function u on S2 (u £ C 2 (S2) n 
C l {S~2)) satisfies the following representation formula 

J u(x)da x - j U*(x,x )A s u(x)da x 
S2 S2 

— j U*(x,xn) curl s u(x) ■ t(x)ds 7 - + J u(x) curl s U*{x, Xn) • t(x)ds 7 - 
c c 

j u(x ) if x Q € Si, , , 

"\0 ifxaeSXSL {Z -*> 

Proof. We obtain Green's second formula by interchanging the roles of u 
and v in (12, ip , adding the two identities and using the symmetry of the left 
hand side. 

u[x)A$v[x) — v(x)Asu(x) da x 

\v{x) curl yitfx) — u{x) curie ; v(x)] ■ t(x)ds x . (2-6) 

We define the e-neighbourhood of xn on S, B e {xp) := {y £ S : \y — xq\ > e} 
and set S2 t£ '■= <S \ B £ (xq). The second Green's formula for S2 <£ with v(x) = 
U*(x, xq) yields: 

u(x) da x - I U*(x, x ) A s u(x) da x 



S-2 



4-7T 

S2,e S2 



\U*(x, xn) cur\g u(x) — u(x) cuvl.q U*(x, xq)] ■ t[x)ds x 

[U*(x, xn) curlg u(x) — u(x) curb ? [7* (x, xp)] ■ t(x)ds x (2.7) 



8B e (x ) 
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First we observe that 



J U*(x,x ) A s u(x)da x \ < ||A 5 «|| L oo( S ) J \U*(x,x )\da x < M, 



S2 S2 
and hence: 



lim / U*(x, xq) Agu(x) da x = / U*(x, xq) A^ufx) da x . 
^0 J J 



S2, e $2 

Furthermore, we can estimate the integral 



J U*(x, xq) curl^ -ufx) ■ t(x) ds x < I [ curl^ nl | t,°° J \U*(x, xq)\ ds x , 

3B £ dB e 

By changing coordinates, without loss of generality xq can be taken to be 
the north pole of the sphere, i.e. xq = (0, 0, 1) T . The curve dB £ is then 



fully described by the latitude 6 £ , say. According to the cosine law we have 

2 

cos 6 £ = 1 — tj- and thus 

/\U*(x,xq)\ ds x = — [ log 1 1 — cos 6 £ I sin 8 £ dip 
47r h-K 



dB e 



1 / £ 2 

= --ey(l- T )logy ^0 (ase^O). 

To analyse the second contribution along dB £ , we again assume xo to be the 
northpole. We then compute 

1 TT*( \ Smd Z 

curLjC/ x,x = — — ^r e v 

4tt(1 — cos# £ ) 

Since the line element on the surface of the sphere is given by 

t(x(ip,6)) ■ ds x ( Vt Q) = d9 eg + s'mOdipe^, 
we deduce that (note the orientation of dB £ ): 

/sin 2 6 £ f° 
u(x)\cuvl s U*(x,xo) • t(x)]ds x = -— — e —- / u{tp,0 e ) dtp, 
47r(l - cos 6^ J J 27T 

dB £ 
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which in the limit as e — > tends to 

u{x ) e 2 (l - t) i x 
iim 5 = u(xq). 

2 

Taking the limit as e — >• in (|2.7p proves the result. ■ 

At this juncture, we draw the reader's attention to the term u(x) da x in 
the representation formula (|2.5p . If it satisfied A^u = in 52, the familiar 
integral representation formula for the Laplacian in 2-D would not involve 
such a term; indeed, the left hand side of the representation formula consists 
of line integrals only if u satisfies the side constraint, f s uda x = 0. This is 
linked to the Gauss constraint. 

3 Layer potentials and boundary integral opera- 
tors 

Having derived an integral representation formula for solutions of the Laplace- 
Beltrami problem in 1S2 in the previous section, we are now in a position to 
reformulate the boundary value problem as an integral equation. 

3.1 Single and double layer potentials 

Following the integral representation derived in Proposition 12.21 we define 
the following two layer potentials: 

• The single layer potential with sufficiently smooth density function 
a: 

(Va)(x) := J U*(x,y)a(y)ds y for x C, 
c 

• and the double layer potential with sufficiently smooth density 
function fj,: 

(Wn){x):= J v(y)[curtsU*(x,y)-t(y)}ds y iorx^C. 
c 

By Proposition 12.21 every solution to the homogeneous Laplace-Beltrami 
equation can be written as the sum of a single and a double layer potential 
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modulo a constant. This is the starting point for the so-called direct bound- 
ary integral approach. However, for the purpose of this paper we follow the 
layer ansatz based on the following observation. 
For x $l C, the single layer potential satisfies: 

A s (Va)(x)=A s J U*(x,y)a(y)ds y = J A s U*(x, y) <r(y) ds y 
c c 

= 4tt / a ^ dSy = ( 3 ' la ) 
c 

under the constraint J a(y)ds y = 0. (3.1b) 
c 

Hence, we may find the general solution of the Dirichlet boundary value 
problem (|1.6p as 

u(x) = / U*(x,y)a(y)ds y + p, (3.2) 



c 



where p € M. is some Lagrange multiplier related to the constraint (I3.1bp . 
Similarly, the double layer potential satisfies the Laplace-Beltami equation 
for x C 

A s (Wfi)(x) = A S J ix{y) [curl^C/* (x,y) ■ t(y)] ds y 
c 

H(y) [A s cmLj U* (x, y) ■ t(y)} ds y 

= J v(y)\cm±s A sU*(x,y)-t(y)]ds y 
e 

= 0, 

without any further constraints on the density fi. We might thus also try to 
look for the solution to (|1 .6j) in the form of a double layer. 



c 



3.2 Jump relations for V and W 

In the previous section, we have only defined the layer potentials for x away 
from the boundary curve. However, in order to align the operators with the 
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given Dirichlet data along C, we need to investigate their behavior in the 
limit as x approaches C. Similarly, if one is interested in solving the Neumann 
problem in which the tangential component of the vectorial surface rotation 
is prescribed along C, one has to investigate the limit features of this quantity 
for the layer potentials. In both cases, there will be certain jump relations 
across the curve C. For the purpose of this paper however, we will restrict 
ourselves to the Dirichlet case. First, consider the single layer potential with 
density a for x C: 



(Va)(x) = I U*(x,x)a(x)ds x 

log[l - (x , x)) a{x) ds x (3.3) 



c 

1 

47T 

c 

The following lemma describes the limit behavior of the single layer poten- 
tial. 

Lemma 3.1 For xq £ C we have: 

(Va)(x ) := lim (Va)(x) = [ U*(x ,y) a(y) ds y 

C 

as a weakly singular line integral and hence (Va) is continuous across C. 

Proof. Fix an arbitrary e > 0. Let xq S C be fixed, and and x € S satisfy 
\x — xq\ < e. Introduce the notation 

C £ ,< ■= {y G C, \y - x \ < e}, C e> := {y E C,\y — x \ > e}. 

Then, if we define 

I £ (x) := / U*(x,y)a(y)ds y - / U*(x ,y) a(y) ds y , 



c 



we can easily show 



I £ = [U*(x,y) - U*(x ,y)]a(y)ds y + U*(x,y) a(y) ds y . (3.4) 
The first integral in (|3.4p vanishes in the limit as x — > xq, i.e. 



lim / [U*(x, y) - U*(x ,y)} a(y) ds y = 0. 
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The second term in (|3.4f) we can bound in terms of the density a: 



Sc. 



U*{x,y)a{y)ds y 



< HU«„(C) / \U*(x,y)\ ds y . 



To finish the proof, note that we can estimate 



U*(x,y)\ ds y < / \U*(x,y)\ ds 



X—tXQ 



\U*(x ,y)\ ds y e -^i 0. 



y eC 
\y-x\< 2e 



y€C 
\y - xo\ < 2e 



Putting these estimates together, we see that lim lim I e {x) = 0, which 

E~ ¥0 X— YX 

proves the assertion. ■ 

The case of the double layer potential is slightly more involved, since the 
limit process reveals a hidden delta function. To see this, consider the double 
layer potential with density fi for x C: 

(Wfx)(x)= I n{x) [curlsU*(x,x)-t{x)] ds$ 



4tt 



A(x, x) 



cos((/j — if) cos 9 sin 9 + sin 9 cos 9 
— sin#sin(<^ — cp) 



t x) 



dSn 



(3.5) 



where x = x(<p, 9) and x = x(ip, 9) and A(x, x) = 1 — cos(y> — ip) sin 9 sin 9 
cos 9 cos 9. 



Lemma 3.2 For xq € C we have: 



(<y^Wri(x ) 



lim (Wfi)(x) 

S2^x^tXQ 



= (#a0(*o) + f 1 -^) m(*o), 

where a(xo) represents the interior (with respect to S2) angle of C at xq. 
For a smooth curve, a = n. The operator (Kfi) is given by the following 
integral expression: 

(K/j,)(x ) = lim (K £ fi)(x ) 

£->0 



lim 

e-s-0 



n(x) [ curie ; U*(x, x) • t(xj] dsx- (3-6) 



i 6 C 
\x — xo\ < e 
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Hence the double layer potential satisfies: 



(3.7) 



where we tacitely assumed the orientation of the tangential vector t along C 
to be in accordance with the orientation of S2 in the sense of Stoke 's theorem. 

Proof. Given e > 0, let x 6 £2 with |[x — a?o|| < e.We introduce the 
notation 

C £) < := {x £ C,\x — xq\ < e}, C £j > := {x £ C,\x — xq\ > e} 

Then, 

(W/i)(x) - (K £ /j)(x ) = I u(x) [ curlg E/*(x, x) - curlg U*(x , x)} ■ t(x) ds* 

n{x) curl g U*(x,x) ■ t(x) dsx 
and the first integral again vanishes as x approaches xq, i.e. 



[curl g U*(x, x) — curl y C/*(xq, x)] • i(x) ds^ 



a;— >xo 



0. 



The second term can be rewritten as follows: 
fi(x) curl y U*(x, x) ■ t(x) ds x 

\p(x) — ju(xo)] curlg U*(x, x) ■ t(x) dsi 



+ M^o) / curlg t7*(a, x) ■ t(x) dg^ 



(3. 



For the first integral on the right hand side of (|3.8p we have the estimate 



\p(x) — fJ-(xo)] curlg U*(x, x) ■ t(x) ds x 



< sup \n(x) — [i(xo)\ / I curl g U* (x, x) • t(x)\ dss 

< M ■ length (C e ) ■ sup |/u(x) — /i(xo)| 
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for some constant M, and hence the integral vanishes in the limit as e — > 0. 
For the second integral in (|3.8p . we define £l e (x) := {x & S2 '■ \x — x\ < e} 
to see 

M x o) / curl y U*(x, x) ■ t(x) dsr. = ttixn) I cml s U*(x,x) ■ t(x) ds x 

dQ e (x ) 



/i(xo) / curl y U*(x, x) ■ t(x) dsn 



\x — Xo[=£ 



Using the representation formula with u = 1 we get 



//( .r (1 ) I 1 - j da £ I - fi(x ) J CWI5 r ? *(.r..r) •/'(.?• )f/.«, ; . 

|x - K | = e 

and without loss of generality we compute the remaining integral with re- 
spect to the northpole to find for all xq: 



lim / curl y U*(x, x) ■ t(x) ds 
e^O J 



ol{xq) 



2vr 



\x — xq\ = e 



Putting the parts together we see that 



lim lim 

e— >-0 x— >xo 



(Wim)(x) - (K e Mx Q j) = (l- ^2) //(zo). 



4 A BIE strategy for solving the Dirichlet problem 

With the single and double layer potentials defined as in the previous section, 
we are now in a position to reformulate the Dirichlet problem for the Laplace- 
Beltrami operator. For the purposes of this paper, we assume sufficient 
smoothness of the data g and the curve C such that all the operators are 
well-defined; these assumptions can be relaxed, and the precise regularity 
and smoothness assumptions necessary are a subject of a forthcoming work. 
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In what follows, however, we assume the curve C is at least C 2 , and that 
the boundary data is smooth. We shall present a numerical example where 
we allow g to be Lipschitz; 

Recall that we wish to find a smooth function u such that 

— A$u = in S2, and u = g onC. (4-1) 

We seek a solution of this equation in terms of a layer ansatz. That is, we 
wish to find a density, a or fj,, so that either 

u := Va or u := W \i (4.2) 

solves 14.11 

Lemma 4.1 If the density a solves the boundary integral equation 

Va = g, on C, and also / ads = 0, (4-3) 

Jc 

then the function u := Va solves ( f^.i[ ) -If the density \i solves 

{ l -I + K){i = g, onC, (4.4) 

then the solution of f^.ip is given by the double layer potential, u := Wfi. 
The proof is immediate from the previous section. 



4.1 An indirect integral equation formulation 

For concreteness, we describe in some detail a variational strategy to solve 
a boundary integral equation, whose solution then can be used to solve 
(|4.ip . We seek a solution u of the Laplace-Beltrami operator with prescribed 
boundary values on C. The solution is assumed to be of the form 

u = Va + p 

where the density satisfies the integral equation 

Va + p = g on C 

along with the constraint 

J ads = 0. 



16 



We can write the weak formulation of this problem in saddle-point form as: 
Find a € % and a multiplier p£l such that 



(Va,x) +p(i,x) = (g,x) 

q(l,a)=0, 



(4.5a) 
(4.5b) 



for any test function x S H and any real constant q. Under the present 
assumptions on smoothness, we set H = C(C), and (•, •) is simply the L 2 — 
inner product along C; we will describe the appropriate Sobolev spaces in 
which to naturally seek a in subsequent work. The angle brackets will then 
represent the duality pairings in L? . 

The discretization strategy is now standard. Let be a partition of C, with 
sub-interval size h > 0. We approximate % by a finite-dimensional space, 
Sh, which is parametrized by the meshsize h; as h 0, the approximation 
error mf Vh€ s h \\u — v^Wh ~> f° r au u £ We then study the discrete 
Galerkin problem: 

Find ah G Sh, Ph £ suc/i i/iat /or ali (x^, }/,) £ Sj, x 1, 



We shall provide an error analysis of this system in a subsequent paper, 
based on the correct choices of Sobolev spaces for the densities a and ap- 
proximation spaces Sh ; at present, we present numerical experiments to 
validate the boundary element strategy. 

4.2 Numerical experiments 

In what follows, we choose C to be the equator of the sphere, which is 
described by the latitude 9 = tt/2. We solve the Laplace-Beltrami equation 
in the southern hemisphere £2- To do this we prescribe Dirichlet data g 011C 
and solve the discrete Galerkin system (|4.6p . The partitions = U^fij of C 
are chosen to consist of uniform sub-intervals of size h. The approximations 
are sought in the space of piece-wise constant functions, i.e. 



In the specific case of the southern hemisphere, a Green's function for the 
problem is known and we can write down the solution in closed form as 
follows: 



{Va h , Xh) + Ph{l,Xh)ds = (g, Xh) 
qh(l,(7h) = 0, 



(4.6a) 
(4.6b) 



Sh ■= {X I X(¥>) = (H for ip € OJ. 



cos(0) 



g(<po)dtpo 



(4.7) 



1 + cos(<^ — ipo) sin(6>) 
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Table 1: L 2 error of BEM solution, measured along 6 = 3/4 it 





Convergence for g = sin(</?) 


Convergence for g = h{ip) 


DoF 


L 2 error 


Ratio 


L 2 error 


Ratio 


20 


2.41e-4 




0.023 




40 


2.99e-5 


8.06 


2.57e-4 


89.49 


80 


3.73e-6 


8.02 


3.06e-5 


8.40 


160 


4.66e-7 


8.00 


3.73e-6 


8.20 


320 


5.83e-8 


7.99 


4.61e-7 


8.09 


640 


7.28e-9 


8.01 


4.73e-8 


8.05 



for 6 > 5 and 92 E [0,27r]. This expression serves as a reference for our 
computed solution. 

We report the convergence behavior of the method in terms of the L 2 error 
of the computed solution along the latitude 6 = 3/4 7r. In Table 1, we report 
this L 2 -error versus the number of unknowns, for two choices of Dirichlet 
data: g = sin(^), and g = h((p), where h is a hat-shaped function, 



h(<p) :- 




otherwise. 



We note that halving the mesh-size reduces the 1? error by a factor of 8 



in both cases. Figure 2(a) shows the convergence behavior in terms of the 
above mentioned error versus the number of unkowns. In figure 2(b), we 



see the actual solution of the boundary value problem, corresponding to the 
piece- wise Dirichlet data g = h((p). 




(a) Convergence behavior for different (b) View on solution from below the 
choices of g south pole for hat-shaped data, h(ip). 
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5 Conclusion 

We have presented a boundary integral formulation, and associated Galerkin 
discretization strategy, for a boundary value problem for the Laplace-Beltrami 
operator on the unit sphere in M 3 . Numerical experiments verify the appli- 
cability of the idea, and a rigorous error analysis will be presented in future 
work. 
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